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1. Introduction 

The equation of state of a quarks gluon plasma at high temperatures and/or 
densities is one of the most important unknowns in our current understanding 
of strong interaction physics 1 . The applications of such an equation of state are 
varied, ranging from cosmological compact objects to the physics of heavy ion col- 
lisions. Unfortunately, due to the high degree of non-linearity present in QCD, 
the determination of this equation of state has proved to be a difficult task. For 
example, even within perturbation theory, infrared singularities require the calcu- 
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lation of an infinite number of graphs for the partition function beyond fifth order . 
General expressions for Green's functions are available for the case where the in- 
ternal momentum is large (~ T) while the external momenta are soft (~ gT), the 
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so-called hard-thermal loop region ; using standard manipulations , one can then 

determine the partition function corresponding to all hard modes in the theory. 

The soft-mode contributions to the partition function have been studied using var- 

ious approximations and numerical calculations have also been developed (though 
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not to the extent as in the T = case) . 

In this paper we propose a new approximation within which the physics of a 
quark-gluon plasma can be studied. The formalism is based on the Thomas-Fermi 
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model of the atom and will be called Thomas- Fermi QCD (TFQCD). We con- 
sider a plasma of quarks and gluons confined to a volume V which we imagine 
subdivided into a number of subvolumes, each of which is large enough for the 
partons they contain to be considered a statistical ensemble. These subvolumes 
interact via a background gauge field whose sources are the thermally-averaged 
non-Abelian charge densities of the subvolumes. The subvolumes are assumed 
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to be small enough for the background field to vary very little inside them, and 
because of this the background field sources are essentially point-like. The require- 
ment of stability, together with the Yang-Mills equations for the background field, 
furnishes a closed set of equations which can be solved; from the solution the equa- 
tion of state for the system can be derived. This program requires the evaluation 
of the thermally-averaged non-Abelian charge densities which we obtain using per- 
turbation theory; in this paper we will use the lowest-order approximation, but a 
systematic improvement is straightforward. 

The atomic Thomas- Fermi approximation is useful when calculating bulk prop- 
erties of an atom with a large number of electrons, such as the total ionization 
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energy ; it is also useful as a starting point for a Hartree-Fock approximation. We 
expect the TFQCD model of a quark-gluon plasma to be reasonably accurate for 
bulk properties of the system, such as the equation of state. There are some differ- 
ences between the atomic Thomas- Fermi and the TFQCD formalisms; in particular 
note that, in contrast to the atomic case, the quark-gluon plasma is not stable: if 
left alone it will fly apart and undergo a phase transition into a gas of hadrons. In 
order to study a gas of quarks and gluons we are forced to imagine the system to 
be enclosed in a container at sufficiently large temperature and/or density. 
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The presence of an external confining agency is reminiscent of the bag model ' . 
Through most of the paper we will consider, in contrast to the usual bag models, a 
situation where the partons are not confined, and for which the external pressure 
is assumed to be generated by a some physical apparatus. Despite this difference 
the bag boundary conditions are also relevant for the present model: the system is 
assumed to be confined to a spherical volume out of which neither fermion number 
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nor color can escape, this requires we impose both the original and chiral bag 
boundary conditions. 

We will also briefly study a system corresponding to a hadron at zero temper- 
ature, and will show that the bag constant and strong coupling constant obtained 
in the present approach are consistent with those obtained using the bag model. 

The volume of the system V will be kept finite in all computations; the results 
will then include finite-volume effects (such as terms in the extensive thermody- 
namic quantities proportional V 2 / 3 ). In the infinite volume limit these surface 
effects can be neglected and the equation of state reduces to that of an ideal gas 
of gluons and quarks. 

In the following section we will describe the construction of the TFQCD model 
and present some simple applications. We will concentrate on the case of an SU (2) 
gauge theory with a single species of massless fermions, and then describe the 
modification required for the important case of an SU (3) gauge theory with three 
(massive) fermion flavors. The remainder of the paper is organized as follows. In 
section 3 we derive the equation of state for this model in the cases of zero baryon 
number and zero temperature. The discussion of the extension to SU(3) and to 
more flavors is presented in section 4. Some parting comments are presented in 
section 5; and a mathematical detail is given in the appendix. 



4 



2. Description of the model 



The model we propose is, as mentioned above, an extension of the Thomas- 
Fermi model to the case of QCD. We consider a gas of partons inside a volume V; 
we then imagine partitioning V into small subvolumes 5V which are big enough so 
that the partons (quarks and gluons) contained in them form a statistical ensemble 
determined by a temperature T and, for the fermions, a chemical potential \i. Each 
subvolume is required to be in equilibrium with its environment which implies that 
the temperature and chemical potential are the same throughout the system (this 
is intuitively obvious, we present a proof in Sect. 2.3). The system is also assumed 
to be static so that no currents are present. 

We assume that the subvolumes have a non-zero average color charge, which 
implies that the zero-component of the gauge field goes to a constant, A at its 
boundary 11 . We will refer to A as the background gauge field. The background 
field is assumed to vary slowly and smoothly between the 5V, and is determined self- 
consistently by requiring it to satisfy the Yang-Mills equations corresponding to the 
average charges of the subvolumes (which themselves depend on the background 
fields). This approach presupposes that the magnitude of the charge in any given 
SV is small, and that the background field is approximately constant within each 
subvolume; both these assumptions will be verified a-posteriori 

Finally, we also assume that our system is spherically symmetric; this re- 
quirement considerably simplifies the calculations yet preserves the essential non- 
Abelian character of the problem. The equations obtained for the backgound fields 
are then similar to the ones derived when considering the coupling of classical, 
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spherically symmetric Yang-Mills fields to external sources ' . 

In the rest of this section we will treat the various ingredients of the model 
separately. We first review the Yang-Mills equations within the spherically sym- 
metric Ansatz. We then obtain the expression for the partonic sources for the 
background fields and the various thermodynamic observables. Next we derive the 
stability conditions for the system. Finally we combine these results in order to 
obtain the equations for the background fields which determine quantitatively the 
Thomas- Fermi-Q CD (TFQCD) model. 

The conventions which we use are the following. The model is based on 
an SU(N) Yang-Mills theory with one species of massless fermion; the (anti- 
hermitian) group generators are denoted by T a and the gauge coupling constant 
by g. The co variant derivative is = + where = gA^T" 1 . The full 
Lagrangian is 

C = i$W>-\(F«,) 2 (2.1) 
where ip denotes the quark field, A the gauge field, and 

F; v = d^A a v - duA% + ge ahc A\A%. (2.2) 

The sources are 

3l = i$T a ^. (2.3) 

Latin indices from the beginning of the alphabet (a, b, c, etc) correspond to color 
indices; Latin indices form the middle of the alphabet (i,j,k, etc) denote space 
indices. 
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2.1. Spherically symmetric gauge potentials and equations 

As mentioned in the previous section we will assume that the long-range forces 
in our system are described by a non-Abelian background gauge field generated by 
the average charge of each subvolume. We also assume the system to have spherical 
symmetry. Thus we need the most general expression for a spherically symmetric 

14 

non-Abelian gauge field, which is well known , and is reviewed for completeness 
below. 

The most general spherically symmetric Ansatz for the gauge potentials of an 
Abelian theory is simply A = (f)(r,t), A = a(r,t) r, where A denotes the vector 
potential and r = |r|. It is clear, however, that we can choose a gauge where 
a(r, t) = 0, so we can take A = 0. 

For the SU(2) non-Abelian case the structure is much richer ; the most 
general spherically symmetric Ansatz is 14 (the over-bar denotes the background 
fields) 

A° a = A r a 

A\ = eiajij ^- — + ($ia - r;r a ) y + r^r^i ^ ^ 

which exhibits spin-isospin mixing . The fields </?i ; 2 and Ao,i depend on r and t. 
Within this Ansatz the SU(2) Yang-Mills Lagrangian becomes 



_pai = _ez \T>$\ + — l$l (2 5) 

where the indices /i, z/, etc. equal 1 (corresponding to r) or (corresponding to t); 



=#=1 The situation is similar for larger groups, see section 4. 

#2 In this respect the present approach differs from other investigations into spherically sym- 
metric hadron physics; see Ref. 15. 



the metric is diag(l, —1). We also defined = d^Av — duA^, $ = <pi + iip2 and 
T^\i — dft + igA^. The above expression is invariant under the gauge transformation 

A* - A* " d^A, $ - e^ A $, (2.6) 

which is a remnant of the original non-Abelian invariance. 

We now consider the coupling of the above fields to a spherically symmetric 
charge density p a , where spherical symmetry requires p a = qf a . The coupling is 
then described by adding a term 

^interaction = fi^oP" = 9 Q A (2.7) 

to the Lagrangian. 
2.2. The partition function for gluons and fermions 

We now imagine that the volume of the system, denoted by V, is subdivided 
into a large number of subvolumes 5V. The gauge fields inside each subvolume are 
separated into a background piece A ^ and a fluctuation a a ^\ A^ = A® + a^. 

In this subsection we evaluate the partition function for the partons inside 5V. 
This object, which we call Z$v, will depend on A^, and we can use this dependence 
to obtain the thermal average of the non-Abelian currents, 

Since the system is supposed to be in a static configuration we require jf = 
which implies we can take Af = inside 5V. Since the background fields are 
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assumed slowly varying, we also take Aq constant inside 5V; we then choose a gauge 
such that Aq is diagonal inside 5V. Hence Z$y will depend on the temperature T, 
the fermionic chemical potential /i, and the n components of Aq associated with 
the diagonal generators (n is the rank of the gauge group). 

As a first approximation we will neglect the interaction between the fermions 
and the a a ^ as well as the non-linear couplings among the a a ^ these interactions 
can be included perturbatively. 

Concerning the scale of 5V we will assume that it is set by the fermion thermal 
wavelength, A, that is, 

5V ~ A 3 . (2.9) 

We will for the moment restrict ourselves to the case where the gauge group is 
SU(2) (the extension to SU(3) will be described in section 4 below). In this case 
the group generators are T a = a a / (2i) and, within 5V, Aq = gA^T® = gAoas/(2i), 
so that 

A b = a b + A 5 b3 , A\ = a\ (inside SV) (2.10) 

We first evaluate the fermionic contribution to the partition function, and then 
calculate the contributions from the a^. 

Fermionic contribution. 

When considering the fermionic partition function we will assume only one 
massless fermions species (the modifications required by several species and/or 
non-zero masses are straightforward) . Thus we look for an approximate expression 
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for Z^ = det (i lp + 7x70) where the gluon fields take the form (2.10). Inside a 
subvolume 5V it is assumed that the fermions behave as a statistical ensemble, 
that the interaction with the is small, and that the background fields A® are 
essentially constant. Adopting these approximations we reduce the calculation to 
evaluating det [i $ + (— iAq + //) 70] with Aq = gAoas/(2i), Ao =constant. 

The partition function for an ideal gas of massless fermions at temperature T 

4 

in a volume 5V and with a chemical potential equal to [i is given by 



lnZ 



(35V 
12^2 



/i 4 + 2(ttA;T)V + t^(^T) 4 
15 



(2.11) 



where the zero subscript indicates that no gauge fields are included. The constant 
background gauge fields are then included by replacing jj, — > /i ± gAo/2 in Zq with 
the sign depending on the isospin of fermion, and adding the contributions from 
each isospin component. Thus, within the above approximations, we obtain 



InZ 



(35V 
6tt 2 



9 A 9 2 Aq 42. ^ rr\A 



(2.12) 



where we defined the thermal wavelength 



A = 



2^2 + (nkry 



(2.13) 



This approximation to Z^ generates the following expressions for the fermionic 
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contribution to the (local) thermodynamic quantities 



q * g(35V V dA ) 



60 2vr 2 [24A 2 48 
1 fdlnZ^ 



2k \g 2 j 



ln^V> _ 1 \g 



(35V 6tt 2 [ 16 



1 \9 2 Al gUt 




(2.14) 



where V denotes the pressure and s, e, n, q, etc denote the entropy, energy, particle 
and charge per unit volume. Note that = 37^ — 2Aoq includes the energy of 
the interaction with the gauge fields. 

Gluonic contribution 

The gluonic contribution to the partition function is obtained in a manner 
similar to the one followed for the fermions. As before we will ignore the self- 
interactions of the fields a^, in this case the partition function, including the 
Fadeev-Popov determinant, reduces to Z g i UO ns = det — D 2 d j where D a dj denote 

16 

the covariant derivative for the background fields in the adjoint representation . 

In calculating this determinant we take into account the partons are supposed 
to be in a box of side ~ A. Moreover, the background field is supposed to include 
the effects from the zero (Fourier) modes in the field. It follows that we need to 
include only modes with energy above 




(2.15) 
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Using the gauge = S^fiS a,s Ao gives 

In Z gluons = -V J ^ [in (l - e -/»(P-**0) + In (l - e "/Wo)) + ln ( x _ e -ft>y 

P>Po 

(2.16) 

which corresponds to a gas of massless bosons with chemical potential ±#.4.0 and 
0. 

We will argue below (section 2.4.1) that the background field Ao is monotonic 
in r and that XgAo < 3ir/2. Using also the fact that (3p > Atc 2 we find that to a 
good approximation 

^giuons = + l^p! (^2 + 2(]po + 2) e -Ppo [sinh {(3gAo/2 )] 2 (2.17) 

This shows that the deviations from the free-gluon values are exponentially sup- 
pressed (recall that (3p > Air 2 ) and can be neglected. In this case 



I giuons — 15^2 «gluons — 5^2 . . 

s gluons — 157,- Qgluons — u 

where P denotes the pressure, and e, s and qthe energy, entropy and charge per 
unit volume respectively. The errors incurred are below a few percent for the 
thermodynamic quantities and below 0.0075% for the charge. 

2.3. Stability conditions 
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The stability criterion can be obtained from the Wong equations , but a 

18 

more elegant argument can be gleaned from a paper by Brown and Weisberger . 
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Consider the background field contribution to the energy momentum tensor #back> 
which satisfies 

^Cck = 9PaK°. (2.19) 

where denotes the field strength for the background fields, and p a the thermally- 
averaged non-Abelian charge density. 

Since the total energy momentum tensor is conserved, it follows that the aver- 
aged partonic contribution 9 pa _ rt satisfies 

^Ct = ~9PaK°. (2.20) 
For static situations the above equation implies 

diO^t = -9PaFi°- (2-21) 

If in addition we impose spherical symmetry (see section 2.1) p a oc r a , Aq oc f a 
which implies p a Fi° = —padjA^. For a homogeneous gas of partons the space 
components of the energy momentum tensor are 6*p art = PpartAj- Collecting these 
results we get djVp^rt = gp a djAQ, or, equivalently 

<tP pta t = 9PadA$ (2.22) 

which is the desired constraint. 

For the Abelian case (2.22) reduces to the usual Thomas-Fermi equilibrium 
condition: the pressure on SV is balanced by the electrostatic force. 
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This stability condition requires the chemical potential and temperature to be 
r independent. Indeed, lnZ^y, the parton partition function for a small volume 
5V, is a function of T, Aq and /i; using (2.14) we obtain 



=gqdA + n P art dfi + (e part - /in part )^ 

(2.23) 

where n pa rt and e pa rt are, respectively, the particle and energy densities of the par- 
tons, and p a = qf a . Substituting this expression for p a , using (2.4), and comparing 
to (2.22) we obtain dp = dT = 0. 

2.4. The TFQCD equations 

The equations of motion are derived from the spherically symmetric Lagrangian 
for the background fields (2.5) when the potentials interact with a source q accord- 
ing (2.7). The source, given in (2.14), is itself a function of the potentials. The 
resulting equations are 

^ , g 2 f^a i 



\ g z J (2.24) 



(r 2 f^) + 2gIm($*V v $) = -gr 2 q5 vfi . 

The gauge invariance of these equations allows us to chose the Ai = gauge. The 
second of the above equations gives, when v — 1 and for static configurations, 
Im$*$' = 0, so that we can choose $ to be purely imaginary. 
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We will use the notation 



A = ZW, $ = - a (r). (2.25) 
rg ig 
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Then, using (2.14), the above equations become 



7y a =0. 



) 



(2.26) 



where a = g 2 /(A-rc) and x = r/X with A defined in (2.13). 

These equations determine the background self-consistently. Their solution 
requires the specifications of the boundary conditions to which we now turn. 

Boundary conditions The conditions near the origin are determined by considering 
the behavior of Wilson loops as r — > 0. We find that singularities arise unless / 
and a 2 — I vanish at r = 0. Using the a <-> —a symmetry (which is a remnant of 
the gauge symmetry) we can then require a — > 1 as r — > 0. The precise manner in 
which / and a — 1 vanish as r approaches zero is determined by requiring that the 
energy should have no divergences at this point. We then obtain 



It is easy to see that the equations of motion (2.26) require / to be concave or 
convex; since we can always exchange / — > — / we can assume that / is concave. 
In particular this implies that / will not vanish for r > 0. From (2.26) it is also 
easy to show that (f/r)' > 0. In contrast a can (and will) have extrema as well as 
zeroes. 

As mentioned in the introduction we assume that the system is enclosed in a 
container which must be spherical due to requirement of spherical symmetry; we 



/, a - 1 = 0(r 2 ) for r -> 0. 



(2.27) 
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denote by R its radius. If the system is to be confined to the region r < R, there 
should be no leakage of fermion number or color into the region r > R. 

The first of these two conditions (fermion number confinement) requires the 
fermions to have zero radial component of the momentum at r = R. This implies 
that in the vicinity of r = R the fermion gas becomes two dimensional. The 
corresponding (surface) charge density a a takes the form 

a a = df a (2.28) 



as mandated by spherical symmetry. Note however that °& does not have a simple 
analytical form, 



8np 2 



Q(u) = 2 J ds In 
o 



[QiPfji + PgAo/2) - Qtfn - (3gA /2)} 

e s + e u- 



2 ir\ , / 12 In 2 
~- i u + — ) tanh I 2 — u 



(2.29) 



where the analytic approximation to Q is accurate to about 0.62%; the derivative 
is accurate to 0.92%. 

In the examples which we consider in detail we will be interested in the limit 
where Aq is large and where (3 — > or fi — > 0. In these cases we have 



a 2 

1 - wJ- 



(2.30) 



We will require the volume charge density in the bulk to smoothly join the 
surface charge density at the surface layer, that is, p a \ = a a at r = R. Thus we 
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impose, qX = d at r = R which, keeping in mind that the solutions produce large 
values of Aq at R, is equivalent to Aq/(4iv) = A„4q/(37t 2 ), or equivalently 

f(R) = y f (2.31) 

It is of course possible to modify this condition by requiring only that, at r = R, 
q\ = s$ for some number s = 0(1), which is equivalent to replacing A — > A/s; our 
results are insensitive to such a replacement. 

To determine the consequences of the second of the above two conditions (color 
confinement) we need the components of the chromo-electric and chromo-magnetic 
fields parallel and perpendicular to r, 

f..E a = — L (£) r a r-B° = — ^ — ~ 1 f a 

ta\i 



9 2 \\xJ ^ <i 2 \ x 2 (2 32) 

1 a ' (c w -caNi 1 /a 



(f x B a ) 1 = -_- eijo f? (f x E a ) 1 = - — ±- €ija f 



g 2 \x g 2 \ 



x 



The first of these relations, together with the previously derived result (f/r)' ^ 0, 

implies that color will leak from the system unless an appropriate modification is 

9 

included. The situation is identical to the one present in the bag model , and the 
solution which we adopt is the same 10 . We will couple our system at the r = R 
boundary to a CP-odd field rj' via a term proportional to the Chern-Simons term; 
this coupling insures that color is confined to the region r < R W . Denoting by F n ' 
the decay constant on the rj', the coupling to this field at r = R are determined by 
the relations 



f ■ E a = — — r • B a 7]'; f x B a = x E a rj' (2.33) 

from which we derive (r • E a ) (r x E a ) + (f ■ B a ) (r x B a ) = 0; in terms of the a 
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and / fields this becomes 



fa(xf - /) + xa\a 2 - 1) = at r = R 



(2.34) 



which is the desired condition. 

Character of the solutions 

The TFQCD potentials / and a are then obtained by solving the equations 
(2.26) subject to the boundary conditions (2.27), (2.31) and (2.34). These so- 
lutions, as well as all thermodynamic variables, will depend on the parameter 
a = g 2 /(4:Ti). In order to specify a we first fix the thermodynamic variables of 
the system, such as the energy and volume; the TFQCD expresses these thermo- 
dynamic variables as functions of a, which is chosen so that the chosen values are 
met. 

When considering the (2.26) we find that, for given values of X and a, there are 
several solutions satisfying the boundary conditions . Of these solutions there is 
a set (we, in fact, found two such solutions) which minimizes fl, the thermodynamic 
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potential at constant pressure and chemical potential , 



where V denotes the total pressure. Numerical studies show that there is no 
cross-over as a changes: each member of the set of solution which minimizes Q 



#3 Concerning (2.34) we know, from the numerical integration of (2.26), that f(xf — /) does 
not vanish, we also find that it is numerically large for the situations we consider in detail. 
It follows that (2.34) can be approximately replaced by the simpler condition a(R) = 0. 

#4 This is reminiscent of the situations found in the case of classical solutions to the Yang-Mills 

13 

equations with external sources ). 




(2.35) 
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is a smooth function of a. Selecting the solution which minimizes Vt we then 
determined a by matching the specified energy and baryon number. 

The explicit expressions for Q, the total energy £ and the total number of 
particles (baryon number) J\f are 

x 

aX 



k 

■— f 2 (2 + f - 

24tt 1 x 2 J I 9tt 



x 



/i 4 + 2fi 2 (nkT) 2 + ^ (nkT) 4 
15 



R A 



;2.36) 



x 



£=- h>-mr- L 



aX 



x 



+ 24tt / l 2 



2 

3tt 



// + 2/i 2 (rckT) 2 + ^ (vrfcT) 4 
15 



(2.37) 



AT 



2/x I 4 

7T | 9 



/j 2 + (nkT) 2 



x 

R 3 + X [ dx f 2 } , 



(2.38) 



where 



X 



R 

A' 



(2.39) 



For future reference we also provide the expression for the (total) entropy of the 
system 



1 2 14,., 
-S = —XirkT I -X 3 A 2 
k 3 13 



x 

+ y fdx 





(2.40) 



Solutions of the equations for / and a can be obtained using standard numerical 
algorithms; due to the singular nature of the equations at the origin the relaxation 
method is best suited. 
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(a) 



Fig. 




(c) 




Examples of solutions f(x) (dashed lines) and a(x) (solid lines) corresponding to 
R = lOfm, £/V = 4 GeV//m 3 and T = 150 MeV (X = 47). Cases (a) and 
(B) corresponds to the solutions which minimize the thermodynamic potential. 
Cases (c) and (d) have larger f2 and represent unstable solutions; the values of a 
corresponding to each solution are indicated. 



(b) 
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Z 1 




1 : 


10 


1 1 — 

20 30 

X 


40 



1 

0.5 


-0.5 
-1 



0.16 



0.08 





z\-0.04 



We present several examples of the solutions in Fig. 1 where we took R = lOfm, 
S/V = 4 GeV//m 3 and T = 150 MeV (which implies X = 47). All the solutions 
in Fig. 1 satisfy the equations (2.26) and the boundary conditions; the solutions 
which minimize the thermodynamic potential corresponds to cases (a) and (b)^ 5 . 

#5 The singular nature of the equations allows for the multiplicity of solutions; we have found 
8 solutions in total (for the given values of £ and V but having different values of a), though 
we cannot assert that this an exhaustive list. Using the relaxation method, the solution 
that minimizes f2 was lest sensitive to the initial trial functions, solutions with larger Q 
become increasingly more difficult to find as the range of initial configurations which relax 
to such solutions of (2.26) becomes more and more restricted. We have not attempted to 
perform an complete study of the properties and number of solutions restricting ourselves 
to finding the one solutions relevant for physical applications together with some examples 
of unstable solutions. 
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These two solutions correspond to indistinguishable thermodynamics (within nu- 
merical errors); for the calculations below this duality presents no complications. 
We have not attempted to study the stability of these solutions against non-radially 
symmetric perturbations 20 . 

Given these results we must now determine whether they are consistent with 
the original assumptions, that is, whether / varies slowly enough to be considered 
constant in a region of width ~ A. We also must determine to what extent are 
color charges screened. The plots presented correspond to both cases (a) and (b) 
in Fig. 1. 



(a) (b) 




10 20 30 40 10 20 30 40 



x x 
Fig. 2 Validity of the Thomas- Fermi approximation, (a): only in the regions near the 

boundary at x = 47 [x > 45.5) and the origin (x < 1) the approximation does breaks 

down, (b): charge is effectively screened throughout the volume (the logarithm is 

base 10). 



The rate of change of / is sufficiently slow provided the potential Ao changes 
little within a region of size A, this is equivalent to 

(//*) ' { } 

a plot of the left hand side of this equation if presented in Fig. 2(a). We see 
that the condition (2.41) is satisfied except in the vicinity of the origin and the 
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r = R boundary. The value of (2.41) near x = presented in Fig. 2(a) is an 
underestimate generated by numerical errors (the equations are singular at x — 0); 
for x - 0, (f/x)'/(f/x) ~ 

The magnitude of the charge in a subvolume <5V ~ A 3 is obtained from (2.14), 
and equals 

^ = A3q =2i^tf + 1 )' <2 - 42) 

A plot of this quantity is presented in Fig. 2(b); as can be seen the magnitude 
of the color charge inside each subvolume is quite small except near the r = R 
boundary: the system does screen its charges quite effectively. 

2.5. Solutions for small X 

When X is small then / will be small also since it is monotonic (this follows 
from the boundary condition (2.31)). In this case the equation for a decouples and 
so does the boundary condition (2.34), 

x 2 a" + (1 - a 2 )a = 0; a(0) = 1, a'(X) [a(X) 2 - l] = (2.43) 

If we define 

a 2 = ±a"(0), (2.44) 

it is easy to see that the solution to the above equation is a function of a,2X 2 . It 
is then enough to assume a 2 = ±1; the general solutions are obtained from these 
by rescaling x. The solutions to the above differential equation (for a 2 = ±1 ) arc 
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presented in Fig 3. The solutions are monotonic, so the boundary condition at 
x = X is satisfied when a(X) = — 1 which occurs only for a<i < 0, numerically 

a 2 X 2 ~ -4.1 (2.45) 

which completely specifies the solution. 
a 
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Fig. 3 Solutions for small X, solid curve: a"(0) = — 1, dashed curve: a"(Q) = 1. 



Again neglecting / and evaluating numerically the integrals gives 

K AO 

£ ~ -— + 4.28(fcT) 4 V 

2°Q1 ^ 
V ~ + i.43(A;T) 4 

a 

For example, at T = 0, £ = 1 GeV, i? = lfm, a ~ 1 and "P ~ 85 MeV/fm 3 . Fixing 
V and T, the coupling strength a drops as 1/5. 

At zero temperature we have J\f ~ 2/i 3 V/(37r 2 ) and 

* " ai2 + 19.74' ' a V 59.04 {lA<) 
so that the equation of state becomes 

a 1.62 
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In this case M = 2, £ = 1 GeV and R = lfm imply a ~ 2.4 and "P ~ 77 MeV/fm 3 . 

The numbers obtained for the case of small X are then quite consistent with 
those obtained using the bag model 9 (except perhaps for a large value for a). Note 
however that in the present calculation the contributions from the non-ideal gas 
terms are very important and the numerical agreement is not trivial. It is also 
true that the present model is far from realistic (being based on an SU (2) gauge 
theory with a single species of massless quarks). These results are therefore quite 
encouraging but not conclusive as to the physical relevance of this model. 

Concerning the other thermodynamic quantities they relapse to their free- 
particle values up to 0(f 2 ) corrections. Note that the adiabats are, in general, 
defined by X =const. which, for the case J\f = 0, imply P 3 V 4 =const. just like a 
relativistic ideal gas. 

When \i = an approximate solution for / which satisfies the boundary con- 
ditions is 

2 

/ ^ 4.57^, (2.49) 

in this case the entropy becomes S ~ l.5kX 3 and the heat capacity equals Cy ~ 
3S; the largest contribution to these quantities (~ 94%) comes from the J f 2 term. 

2.6. Solutions for large X 

In order to study the solutions to (2.26) for x finite but X — > oo it proves 
convenient to define y = x/X. We are then interested in the small y behavior of 
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the solutions and a power series is appropriate, 



f = f2 



2 I 2a 2 \ 4 , ( ^4 /I 5m/| 2wa 2 14m 2 » 6 



(2.50) 

where w = aX 2 /(607r). Numerical simulations indicate that neither / 2 nor a 2 
increase with X which, using (2.31), (2.34), leads to / ~ 37nr 2 /2X 2 and a ~ 
1 — x 2 /X 2 for Thus, as X — > oo, / — > and a — > 1 for re finite. 

For x<X the boundary conditions require a ~ and / 3> 1; the first of the 
equations (2.26) can then be approximated by 



/"-^; (2-51) 



whose solution (using (2.31)) reads 



/ ^ ; (2.52) 

4+ (X-^v 7 ^™ 



Using these results we can evaluate the various thermodynamic quantities for large 
X. For example, 



£-£idealgas ^ ( 2 -53) 

from which we find \£ — £ideal gasl/^ideal gas ~ 1/X. All other thermodynamic 
quantities exhibit this behavior: for large T and fixed R (corresponding to large 
X) the system approaches a mixture of ideal gases. 
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We emphasize that this is not a result of asymptotic freedom (when the running 
of the coupling is included the large X behavior will acquire logarithmic correc- 
tions), but a property of the solutions to the differential equations. In the infinite 
volume limit the charges are screened which requires Aq = (see (2.42)). 

It is also worth noticing that (2.53) explicitly displays the finite-volume cor- 
rections to the ideal gas results. 

3. Applications 

We now consider some applications of the above formalism. We first study a 
system with vanishing baryon number (corresponding to fi = 0), and then consider 
the case of zero temperature. 

3.1. Zero baryon number 

This situation is believed to be of relevance in relativistic heavy ion collisions, 

21 

such as those to be produced at RHIC , where, in the standard picture, the nuclei 
will go through one another leaving behind a region of hot quark-gluon plasma with 

22 

zero baryon number 

The requirement Af = in (2.38) corresponds to ji = which simplifies some 
of the expressions. In particular the only scales in the system are the temperature 
and the volume. The plot of the equation of state is given in Figs. 4 and 5. 
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Fig. 4 Equation of state within the Thomas-Fermi approximation for the Af = case. 
The graph displays the pressure as a function of the volume for several values of 
the temperature (P in MeV 4 , V in fm 3 , T in MeV; the logarithms are base 10).. 
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Fig. 5 Three-dimensional rendition of the equation of state for Af = 0; (P in MeV 4 , V 
in fm 3 , T in MeV; the logarithms are base 10). 



We have determined a by requiring the solution to minimize the thermody- 
namic potential Q when the energy density equals 4 GeV/ fm 3 at T = 150 MeV, R 
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lOfm (which is consistent with the expectations for RHIC); in this case a ~ 1.568. 

If we now allow the system to expand adiabatically, we can use the above 
expressions to obtain the relationship between T and R corresponding to this 
process. This isentropic transformation describes (in an admittedly oversimplified 
manner) the expansion of a quark-gluon plasma. The entropy is gotten from (2.40) 
by setting /i — 0, the result is 

x 

= — X 3 + - f dx fix). (3.1) 
k 135 3J K J K J 

o 

Since S is a function of X only (a consequence of having only two scales in the 
problem, R and T), the equation for the adiabats is X =constant, or, equivalently 
VT 3 =const. corresponding to an adiabatic index 7 = 4. Note that the J f 2 term 
in S modifies the usual free fermion gas relation S oc T 3 ; the corrections are ~ 20% 
(which is smaller than the corresponding contributions in the case of small X, see 
section 2.5). 

We can also easily determine the energy density for this isentropic process. 
From the expression for the total energy in (2.37) it follows that S\ is a function 
of X only (for the /i — case). It follows that at constant entropy £ scales as 
T. The energy density then will scale as T/i? 3 oc T 4 , just as for an ideal gas of 
massless particles. 

Using the expression for S we obtain the heat capacity at constant volume, 

In calculating this expression one must remember that the boundary conditions to 
the TFQCD equations depend on X, so that we should in fact write / = f(X; x); 
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when the partial derivative is taken in (3.2), / must also be differentiated under 
the integral sign. 

3.2. Zero temperature 

We now turn to the case of zero temperature; the dimensional quantities in the 
system are now \x and R. In this case all dimensionless quantities such as S/li will 
be functions of Rli only. The chemical potential is determined in terms of R and 
Af using (2.38) but this must be done numerically since the non-ideal gas term is 
significant and cannot be ignored. The plot of the equation of state for this case 
is given in Fig. 6. 




The equivalent contour plot for various values of N is presented in Fig. 7 
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The equation of state (for the range of variables presented in figure 7) is well 
represented by the equation 

PV 4/3 = z(N); log z(N) ~ 12.82 + 5.46 

The PV 4/3 behavior is a result of simple scaling arguments and is therefore present 
here as well as for small X. In contrast, the M dependence of the equation of state 
is radically different (cf. (2.48)). 

We have also determined the chemical potential as a function of temperature 
and volume. The result is presented in Fig. 8. 
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log 



1.91 x 10 4 



(3.3) 




As T — ► the entropy goes to zero linearly, 



X = 2V3Rfi, (3.4) 



since the fermionic contribution dominates in this limit; we then also have Cy = S. 



S 



nkT 



V3 
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4. Extensions of the method 



The inclusion of more flavors is quite straightforward, the charges generated by 
each simply add. Possible computational difficulties arise when the fermion mass 
cannot be neglected (as is the case for the strange quark) for in this case a closed 
form for the fermionic partition function is not available. We will not pursue here 
this situation further as it involves no new concepts. 

A more interesting extension is obtained by considering SU(3) as the gauge 
group. In this case there are two important modifications. First, within each 
subvolume 5V, though we still have A$ =constant and diagonal, this now implies 
Aq = g(Ao\3+Bo\8) I (2i) . In general Bo ^ 0, so in this case we will have additional 
contributions depending on this new potential. The TFQCD equations are derived 
in the same way as for the SU (2) case. Therefore the presence of the gauge field 
is summarized by the replacements 




(4.1) 



/i - gB /V3, 



in Zq (Eq. (2.11)) The resulting fermionic partition function is then 



In Z, 



* 4tt 2 24 





(4.2) 



+ 
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Using this result we obtain the charge densities 



<23 



<2s 



gAp 

"6tt 2 
9 



£ {Al + B 2 ) + V3 gfiBo + -1 



4A 2 



V3 



Bo + Bo 2 ) + \g» {Al - Bo 2 ) + ^ 



(4.3) 



The second modification concerns the form of the spherically symmetric Ansatz 
for the background gauge potentials. For SU (3) a possible Ansatz takes the form 
(now including a contribution in the As direction) 



A° a = A r a ; 



A a — tiajTj 



A° 8 = B . 



(a = 1,2, 3) 

± + + (5 ta - r,r a ) ?± + fjfa^i; (a = 1, 2, 3) 
r I r 



(4.4) 



Note however that the choice of the SU{2) subgroup in which the potentials 
Aa, (a = 1, 2, 3) reside is arbitrary, and that it costs no energy to change from one 
such subgroup to another; these degrees of freedom are included through a set of 
collective coordinates . The full Ansatz we use is then (we define = \ n A^/{2i) 
where the X n denote the usual Gell-Mann matrices) 



A» -> A» = U^A^U 



(4.5) 



where the An are given in (4.4) and U is a time dependent SU (3) matrix. 



The Lagrangian for the background gauge fields then becomes 



AtrF 2 

2 11 \xv 



\ tr F 2 V + 2 tr F 0t [A h JR] + tr { [A*, JR] [A\ JR] } 



(4.6) 
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where 



IR= UU\ 



(4.7) 



and F is the field strength corresponding to A. When the form of the gauge poten- 
tials in the SU(2) subgroup takes the form (4.4), IR should have no components 
along the generators of the SU{2) subgroup generated by Ai,2,3, that is we take 



8 1 

M = E 2i XnJR: 



(4.8) 



n—4 



: #6 



which considerably simplifies (4.6). The corresponding action is 

R 

S = J d A x^ tr F 2 „ + ±c 2 J dt tr l/W; c 2 = -jj- J dr (a - l) 2 . 



(4.9) 



Numerically the coefficient c can be very large (for the numerical solutions pre- 
sented c ~ 3.5 x 10 3 A) 



We will use the notation 



Bo 



h(r) 
rg 



(4.10) 



whence the TFQCD equations become (a prime denotes a derivative with respect 
to x — r/X), 

'f + h 2 ( . r: , \ h 



X z 07T 

a 2 -/ 2 -! 



+ (aVs v) - + 1 



(4.H) 



x' 



h"=^- 



a 

6/T 



h(f + h 2 ) 



+ 



f 2 -h 2 



+ h 



x 



#6 The simplicity of this result is a consequence of the fact that U is made to reside in 
SU(3)/ SU(2) ~ S 5 , a five-dimensional sphere, where the number of invariants is very lim- 
ited. The solutions of to the classical equations of motion for U are geodesies representing 
a motion along the great circles of S 5 . 
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which can be solved using the same methods as before. Note that h = is not 
allowed when /i ^ 0. 

For the interesting case /x = 0, h = is a solution to the above equations. 
Hence, for zero baryon number, the previous solutions also satisfy the SU (3) 
TFQCD equations. It does not follow, however, that these solutions again mini- 
mize the thermodynamic potential. Note also that even in the case h = there 
is an additional contribution to the thermodynamic functions from the collective 
variables U 

We will not pursue this case further in this paper. A realistic investigation of 
the SU (3) case requires we include mass term for the (strange) quarks, and also the 
contributions of the collective coordinates to the thermodynamics of the system. 
We will consider these issues in a forthcoming publication. 

5. Conclusions 

We have presented an approximate treatment of QCD based on the same ideas 
as the Thomas Fermi atom. Within this framework the thermodynamics of the 
system can be derived and the results can be compared with the experimental 
results which will soon be available. 

The method is based on a subdivision of the system into subvolumes which are 
still large enough to be considered statistical systems. These subvolumes interact 
through an average gauge field whose sources are the thermodynamically averaged 
non-Abelian charges for the subvolumes. These charges, though small, are not 
completely screened due to the assumed smallness of the subvolumes. 
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The formalism was developed in this paper for the simplified case of an SU (2) 
gauge group, though we did provide a brief discussion of the modifications required 
for and £77(3) theory. We also ignored fermion masses and all interactions between 
the partons inside each of the subvolumes. Nonetheless we found that the numerical 
values for the pressure in the small Af case are in rough agreement with the bag- 
model calculations. 

For large temperatures, or densities (I » 1) the solutions to the equations of 
motion are such that all thermodynamic quantities approach those of a mixture of 
ideal gases, with 1/X measuring the deviation from this limiting behavior. This 
feature is not related to asymptotic freedom but a result of screening. 

In the limit R — > oo we have / = and a — 1, and the equation of state 
reduces to that of an ideal gas. This model then provides an approximation to the 
finite-volume corrections to the ideal gas, this is explicitly demonstrated in (2.53) 
which gives the surface corrections to the energy of the system. 

A realistic calculations must be performed for an SU(3) gauge theory with 
massive fermions; the partition function inside each subvolume should be evaluated 
to the highest order available (or possible) in perturbation theory. The inclusion 
of radiative corrections will induce, among other things, a dependence of the (now 
running) coupling constants on the temperature and chemical potential. For the 
present calculation no such effects were included. Finally one should also include 
finite volume effects as well as the corrections induced by the gluonic partition 
function. We will investigate such realistic situations in a forthcoming publication. 

We found two solutions to the equations of motion satisfying the boundary 
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conditions and which minimize the thermodynamic potential f2. Both lead to the 
same thermodynamics and appear indistinguishable except near the origin (at least 
within numerical errors). A complete study of the behavior of these solutions under 
non-spherical perturbations along the lines of Ref 20 is required to determine the 
one which is most stable. We have not performed such an investigation since the 
presence of two such solutions does not alter the thermodynamics derived within 
the TFQCD approach. 

The above treatment was not based on a semi-classical expansion of the par- 
tition function for the complete system. It is indeed possible to consider such 
an approach and use (2.12) as an approximation to the fermionic contribution. 
Then the integration over the gauge fields can be approximated by a saddle point 
method. We have not done this because the effective action which is to be mini- 
mized in the last step is, due to the Thomas-Fermi approximations used to obtain 
Z^/j, unbounded from below. It is found that the solutions will soon violate the 
Thomas- Fermi conditions and the method is not consistent; this is displayed ex- 
plicitly in the appendix for the case of QED. In contrast, the approach described 
in the above is consistent with the original approximations. 



In this appendix we present a semi-classical calculation of the partition func- 
tion of QED using the Thomas-Fermi approximation for the fermionic partition 
function. The general expression is 



APPENDIX 




(A.l) 
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where S = S g + S^p, the first term denoting the gauge contribution, the second all 
terms involving the fermions. By definition we have 



Zip 



J [dij}][dip]e s * (A.2) 
which is approximated by ~ J d^xV^, where 

(A.3) 



YlTX 1 



(pi + e0) 4 + 2 (nkT) 2 (p: + e0) 2 + J- (irkT) 4 
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where denotes the electrostatic potential, and e the charge of the fermions (only 
one flavor is considered). Note that is positive definite. 

Assuming spherical symmetry the gauge potentials are of the form, = 0(r), A - 
a(r)v. Choosing the a = gauge gives the following expression for the partition 
function Z = f[d(f>] exp^g, where 

f3 R 

S eS = 4njdt J dr r 2 [-± (<t>'f + (A.4) 



where j3 denotes the inverse temperature, R is the radius of the spherical vessel 
containing the system, a dash denotes a derivative with respect to r, the radial 
coordinate, and t denotes the Euclidean time variable. 

The integrand in S e g is not positive definite. Consider for example = 
00 cos(kr + v) for constant v. If 0o is sufficiently small and k sufficiently large, 
the first term in S e g will dominate; the larger k, the more negative S e g becomes. 
The problem in this case is that these expressions for the scalar potential violate 
the Thomas- Fermi condition which requires A0'/0 1. this shows that a semi- 
classical treatment of the partition function is inconsistent with the Thomas-Fermi 

38 



approximation. We have verified that the same problems arise in the non-Abelian 
case. 
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FIGURE CAPTIONS 

1) Examples of solutions f(x) (dashed lines) and a(x) (solid lines) corresponding 
to R= lOfm, S/V = 4 GeV//m 3 and T = 150 MeV (X = 47). Cases (a) 
and (B) corresponds to the solutions which minimize the thermodynamic 
potential. Cases (c) and (d) have larger Q and represent unstable solutions; 
the values of a corresponding to each solution are indicated. 

2) Validity of the Thomas-Fermi approximation. Only in the regions near the 
boundary at x = 47 (x > 45.5) and the origin (x < 1) the approximation 
does breaks down. 

3) Solutions for small X, solid curve: a 2 = —1, dashed curve: 02 = 1. 
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4) Equation of state within the Thomas-Fermi approximation for the M = 
case. The graph displays the pressure as a function of the volume for several 
values of the temperature (P in MeV 4 , V in fm 3 , T in MeV; the logarithms 
are base 10).. 

5) Three-dimensional rendition of the equation of state for J\f = (P in MeV 4 , 
V in fm 3 , T in MeV; the logarithms are base 10). 

6) The equation of state within the Thomas-Fermi approximation for the case 
zero temperature case; (P in MeV 4 , V in fm 3 ; the logarithms are base 10). 

7) The equation of state within the Thomas-Fermi approximation for the case 
zero temperature case; (P in MeV 4 , V in fm 3 ; the logarithms are base 10). 

8) The chemical potential as a function of volume and baryon number at zero 
temperature; (P in MeV 4 , V in fm 3 , \x in MeV; the logarithms are base 
10). 
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